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Numerical codes based on a direct implementation of the standard ADM formulation of Einstein's 
equations have generally failed to provide long-term stable and convergent evolutions of black hole 
spacetimes when excision is used to remove the singularities. We show that, for the case of a 
single black hole in spherical symmetry, it is possible to circumvent these problems by adding to 
the evolution equations terms involving the constraints, thus adjusting the standard ADM system. 
We investigate the effect that the choice of the lapse and shift has on the stability properties of 
numerical simulations and thus on the form of the added constraint term. To facilitate this task, 
we introduce the concept of quasi well-posedness, a version of well-posedness suitable for ADM-like 
' systems involving second-order spatial derivatives. 

o 

in ■ I. INTRODUCTION 

In numerical relativity, Einstein's equations are usually solved as an initial value problem. That is, the spacetime 
is foliated with spacelike hypersurfaces. These hypersurfaces, or slices, are characterized by their intrinsic geometry 
^sj (spatial metric gij) and extrinsic curvature Kij. Subsequent slices in a given foliation are connected via the lapse 
function a and shift vector /3* Under this framework, the components of Einstein's equation naturally separate into 
' constraint and evolution equations for the dynamical variables gij and Kij. Thus, a typical procedure to construct 
] a spacetime consists of first specifying Cauchy data {gij, Kij) that satisfy the constraints on the initial slice and 
> then applying the evolution equations to update these data into the next slice. This procedure is known as free or 
unconstrained evolution. 

. In order to propagate {gij, Kij) using the evolution equations, one must provide in addition a prescription for choos- 
' ing the lapse function and shift vector. One generally regards the lapse-shift prescription as a choice of coordinates. 
T-H At the continuum level, given a lapse-shift choice and a specific set of Cauchy data, the evolution equations yield not 
only a unique spacetime metric expressed in a given coordinate system but also evolved data {gij , Kij ) that continues 
Q to satisfy the constraints. 

O^' Numerical approximations are likely to complicate the picture described above. For instance, because of truncation 
jjl^ I errors, the data {gij, Kij) do not satisfy the continuum Einstein constraints but rather their discrete approximations. 
bJ[). Even if consistency between the numerical approximations of the evolution and constraint equations is achieved, the 
evolved data {gij, Kij) will at best satisfy the constraints up to truncation errors. Of course, this is what one should 
expect. However, the presence of numerical errors could also trigger fast growing modes that render the numerical 
' evolution unstable. Characterizing and controlling these growing modes has been and continues to be one of the most 
^ , difficult and demanding tasks in numerical relativity. 

- - 1 The origin of these destabilizing modes is multiple. Unstable modes could be purely numerical (e.g., the Courant 
instability) or they could be present already at the continuum level. An example of continuum instabilities are the 
so called constraint- violating modes Q . These are solutions of the Einstein evolution equations that do not satisfy 
the Hamiltonian and momentum constraints. The general perception is that this class of solutions consists mostly 
of rapidly growing solutions, although there is no formal proof that this is indeed the case. In summary, even if one 
were to start an evolution with "perfect" initial data that satisfies the constraints, truncation errors could cause the 
numerical solution to drift into a constraint-violating, and perhaps unbound, solution. 

Various groups have accumulated numerical evidence supporting as a possible cause of instabilities the particular 
3+1 form used to recast Einstein's equation as an initial value problem. These observations have partially motivated 
the development of hyperbolic formulations of Einstein equations ■ Hyperbolic formulations have the advantage 
that mathematical tools exist to prove existence and uniqueness of the solutions as well as well-posedness, thus 
providing important information for the implementation of stable discretization algorithms. If excision is used to 
handle the singularities, numerical simulations of black hole spacetimes require the imposition not only of outgoing, 
radiative boundary conditions far away from the holes but also of conditions on the excision boundary that respect 
the propagation of physical modes to flow down into the hole. Because hyperbolic formulations also yield information 
regarding the characteristics of all the field variables in the system, in principle this knowledge could be extremely 
useful when applying these boundary conditions. 



1 



Unfortunately, so far, none of the numerical efforts based on hyperbolic formulations of Einstein equations have 
been able to demonstrate their clear superiority. Currently, the formulation originally developed by Shibata and 
Nakamura [|| , and later re- introduced by Baugarte and Shapiro , (BSSN) seems to be the least prone to instabilities. 
Interestingly enough, the BSSN formulation is not explicitly hyperbolic. Just recently, a numerical implementation 
of the Einstein-Christoffcl hyperbolic system ||^ has been applied to single black hole spacetime and produced three- 
dimensional evolutions with stability properties comparable to those using the BSSN formulation 11 . The common 
lore these days is, however, that the standard Arnowitt-Deser-Misner (ADM) ||l^ formulation is the one which most 
easily suffers from instabilities |13-1^]. On the other hand, the ADM formulation has the advantage of containing a 
minimal number of equations, making it attractive for numerical studies. 

The main objective of this paper is to show that, given a choice of lapse function and shift vector, it is possible 
to obtain long-term stable and convergent simulations using the standard ADM formulation if "appropriate" terms 
involving the constraints are added to the evolution equations. Of course, here the key ingredient is to understand what 
constitutes an appropriate choice of constraint terms. We show that at least in spherical symmetry, it is possible to 
obtain a definite prescription for adding these constraint-terms that only depends on the choice of gauge or coordinate 
conditions, namely the lapse and shift. The idea of adding constraint-terms to the evolution equations (adjusted ADM 
formulations) is not new [|l6|. Just recently, Yoneda and Shinkai |l^ studied in full detail the propagation of the 
constraints in the family of adjusted ADM systems. 

In order to help identify what constitutes a suitable choice of constraint-terms, we introduce a definition of quasi 
well-posedness. The idea of quasi well-posedness is simply to make choices of lapse and shift such that, when the 
ADM system is enlarged with the addition of new variables (typically, but not always, quantities related to the first 
spatial derivatives of the three metric), one obtains a strongly hyperbolic system. Once a strongly hyperbolic system 
is obtained, one can make use of the standard results that apply to those systems regarding the existence, uniqueness 
and well-posedness of the solutions. The essence of quasi well-posedness is then that the properties of the strongly 
hyperbolic system are to some degree inherited by the original ADM system. The advantage of this procedure is that 
by looking at the strongly hyperbolic system, one is also able to gain insight about the propagation of the constraints 
and characteristics of the original system, thus improving the chances of obtaining stable numerical evolutions. 

In Sec. we introduce the concept of quasi well-posedness. The explicit form of the equations and initial data 
used in our numerical evolutions are described in Sec. III. The prescription to adjust the standard ADM equations is 
presented in Sec. [V. We investigate in Sec. ^ the implications of the presence of gauge modes. Next, in Sec. VI, 



discuss the quasi well-posedness properties of the adjusted ADM system under three lapse and shift conditions. For 
each of t hese conditions, we present numerical results showing the stability properties of the simulations. We end in 
Sec. VII with concluding remarks. 



II. QUASI WELL-POSEDNESS 



To understand the general idea of quasi well-posedness, let us consider the initial value problem for the wave 
equation du V = dxx ignoring boundaries. By adding a new variable, n = dt ^p, this equation takes the form 



dtU = P dx 



(2.1) 



where 



P 





1 



and Q 



1 




One can view the system of equations (2.1) as the analogue of the standard ADM equations. Now, adding one more 
variable, = dxf, we rewrite (|2.lD as a first order system in both time and space. That is, 



with 





dtu -- 


^Adx 




-Bu , 






f 


0^ 













and B 






lO 1 







(2.2) 




Notice that one needs to supplement the system (2.2) with the the constraint C — — dx ip- For this simple case, one 
can find estimates for the solution u fro m (2.1) directly, but the idea here is not to do so but rather to use what is 
known about the first order system (2.2). The matrix A has the following set of eigenvalues and eigenvectors: 
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Ai = 1 with eigenvector ei = [0, 1, 1] 

A2 = with eigenvector 62 = [1,0,0] (2.3) 
A3 = — 1 with eigenvector 63 — [0, —1, 1] . 



Given that the matrix A has real eigenvalues and a complete system of eigenvectors, the system of equations ( |2.2| ) 
is strongly hyperbolic |^ and thus well-posed. In this case, the matrix A is also Hermitian (symmetric hyperbolicity) . 
and has real and dist inct eigenvalues (strictly hyperbolic) 



Since the system ( |2.2| ) is well-posed, we can establish existence and uniqueness of the solution to (p.2[), together 
with the usual condition: 

Mt, )|p<De'^*|KO, )|p, (2.4) 

with a and D independent of the initial data and || || the L2 norm. Furthermore, because the vector u contains u as 
part of its components, 

Mt, w<mt, )|p. 

It follows then from equation ( |2.4| ) that the solution of the system ( |2.1[ ) satisfies the condition 

\\u{t, )|p <i^e-*(|KO, W + \\d,u{0, )||2) . (2.5) 



We should remark that, although (2.5) was derived using L2 norms, the solution could in general be bounded by 



the initial data in some other Sobolev norm. The importance of the condition (2.£) and similar bounds is that 
it guarantees, among other things, that the solutions to (^|^) will not have unbounded high frequency modes. In 
numerical calculations however, this is not enough. For practical purposes, even large power-law growths such as 
are likely to be extremely problematic |^ . 

The idea is then to investigate quasi well-posedness in the standard ADM formulation. That is, given a choice 
of lapse and shift, one enlarges the standard ADM system introducing new variables that render the system at 
least strongly hyperbolic. Existence, uniqueness and well-posedness for the ADM equations will then follow from 
the corresponding usual properties of the strongly hyperbolic system. In some instances, to obtain such hyperbolic 
formulations it is not enough just to add the spatial derivatives of the metric as new variables, but one also has to add 
the constraints to the evolution equations, or even take combinations of these evolution equations. But the general 
idea is the same: add the constraints and/or take the same combination in the ADM equations as one does in the 
strongly hyperbolic system. It will be in this sense that we will refer, in what follows, to the quasi well-posedness of 
some ADM systems. 

It is important to mention that for this analysis we have neglected the effects from boundary conditions. For the case 
under consideration (single black hole evolutions via excision of the singularity) boundary conditions at the excised 
region do not affect quasi well-posedness since, as we shall see, all the characteristics of the system are outgoing into 
the hole. At the outer edge of the computational domain, we use the exact analytic solution as boundary condition. 
Our numerical results show that this does not influence the properties drawn from quasi well-posedness for the class 
of gauge choices and added constraint-terms we considered. 



III. THE G - k EQUATIONS 



Although the idea of quasi well-posedness can be applied to more general systems, throughout this paper we will 
concentrate on the ADM evolution equations {g — K equations) obtained from the vanishing of the Ricci tensor, i.e. 

dtgij - Cpgij = ~2a Kij , (3.1) 
dtK,, - = -V,; V,a + a (i?,, + K K,, ~ 2K,k K'' j) (3.2) 

As customary in numerical relativity, we will focus on free evolutions; that is, we will not enforce the constraints 



^Following |18|, a system is called symmetric hyperbolic if the matrix A is Hermitian. It is called strictly hyperbolic if the 
eigenvalues are real and distinct; it is called strongly hyperbolic if the eigenvalues are real and there exists a complete system 
of eigenvectors; and, finally, it is called weakly hyperbolic if the eigenvalues are real. 
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2H = R + K"^ - K,jK^^ =0 



(3.3) 
(3.4) 



except at the initial data. Finally, to facilitate our analysis, we will only consider a single black hole in spherical 
coordinates. The most general time-dependent spherically symmetric metric in a 3+1 form is: 



ds^ = (-a^ +0^/3^) dt^ + 2a^ Pdtdr + a^ dr^ + dn"^ , 



where all functions are assumed to depend only on r and t; dfl^ 



d0^ 



sm 



coordinates, the spatial metric, extrinsic curvature and Ricci tensor are diagonal: 

gtj = diag(a^, 6^, 6^ sin^ 6*) 
X% =diag(7^„ Kb, Kb) 



The g-K equations take the form 

{dt - Pdr) a 
(dt-l3dr)b- 

{dt-Pdr) Ka 
{dt-Pdr) Kb 

and the constraints 



j = diag(i?a, Rb-, Rb) ■ 



-a a Ka + adr P 
-ah Kb 



— drCL drCt 

a 



2 drbdra + a [Rb 



1 

ba 



[Ra + {Ka + 2Kb) Ka] 



{Ka + 2Kb) Kb 



M = drKb 



Rb + 2KaKb + Ki = 

+ iKb-Ka)^ = 0, 
b 



where the components of the Ricci tensor are given by 

2 



Ra 



Rb = 



aH 
1 



a 



3 1,2 



-ad^b + drtt drb) 
—bad^b + bdrttdrb + — a {drb)" 



(3.5) 

and j3 = . In these spherical 

(3.6) 
(3.7) 
(3.8) 

(3.9) 
(3.10) 

(3.11) 
(3.12) 

(3.13) 
(3.14) 

(3.15) 
(3.16) 



Our numerical results consist of free-evolutions of analytic, single black hole initial data (a, 6, Ka, Kb) in a com- 
putational domain < r < Vo, where Tq denotes the location of the outer boundary far from the hole and Ve the 
excision boundary inside the black hole horizon. The numerical code used solves the g — K equations in the interior 
of the computational domain by the method of lines. That is, the g — K system has the generic form 



dt u = P dr u + p , 



(3.17) 



with u = u{a, b, Ka, Kb) and p given by the r.h.s. of equations (3.9-3.12). The starting point is then to approximate 
the derivatives appearing in p by second order accurate, centered finite differences. On the other hand, the the spatial 
derivative in the "advection" term /? drU is approximated by an upwind discretization pi| of the form: 



(Mj 



-1-1 



2 A 



iUj -t- 3Ui±l ~ Ui±2) 

3A^^ 



(3.18) 



where g > and the choice of sign is given by ± = sign(/?). As we shall later see, for th e blac k hole spacetime metrics 
under consideration, the shift vector is always non-negative. Thus, the discretization (3.18) is only needed with the 
upper sign. The above upwind discretization has truncation error: 



dr 



1 



= ^ (1 - 2 q) Ar' diu,Tl df- 



6 



6 



(3.19) 



Thus, the discretization is second order accurate for any choice of the adjustable parameter g, except for q = 0.5 when 
the truncation errors are of 0{Ar^). In our simulations, we set q = 0.5. This value yields both stable evolutions and 
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minimizes the a moun t of dissipation introduced by the discretization. After the discretization of the spatial differential 



operators, Eq. ( 3.17 ) becomes a set of ordinary differential equations for the interior grid points {ri}i=i^,,,^N-i- The 
grid values tq and r^r denote the locations of the inner and outer boundaries, respectively. The temporal updating of 
these equations is carried out via an iterated Crank Nicholson method pOl . Notice that at the outer- most grid point 



tat-i, the upwind scheme ( 3.18 ) requires "ghost" values at r^v+i. These values are computed using second or third 
order accurate extrapolations. 

Next is the implementation of boundary conditions. As mentioned before, there are two boundaries in the problem, 
one far from the black hole a,t = ro and a second at the point of excision, tq = rg. The boundary condition 
used at rg is that all field variables take the values provided by the exact analytic solution. It is important to notice 
that in previous efforts in numerical evolutions of spherically symmetric black hole spacetimes, the stability of the 
evolutions depended not only in the location of the outer boundary but also on the imposition of "outgoing" boundary 
conditions. One of these outgoing conditions is constructed by assuming that the solution error u — u is of the form 
e{t — r)/r" , with u the analytic exact solution. This condition yields 

dtu = -dr u + dru + n ^ !- , (3.20) 



to be applied at the nodal grid point — Vo The outgoing condition (3.20) is then added to the system (3.17) and 
handled with the same method of lines used for the interior points. Obviously, there is no reason to believe that the 
solution error should behave as e(t — r)/r" for general choices of lapse and shift. Another approach to construct outer 
boundary conditions that has been used is to blend [ pl| the numerical solution to the analytic one beyond a certain 
radius, r. 

Finally, at the inner or excision boundary, the working assumption is that all of the fields have outgoing (into the 
hole) characteristics. We will show that for the gauge or coordinate conditions under consideration, this is indeed 
the situation. Therefore, there is no need to impose any boundary condition, and one can just apply at Tg the same 
system of equations used in the interior of the computational domain. However, finite difference discretization would 
require "ghost" values at r_i. We construct these values by extrapolation. Another possibility, which is becoming 



popular in three-dimensional simulations ||22| , is instead to extrapolate the r.h.s. of Eq. (3.17) to ro ~ r^. Either 
approach was stable for the cases considered in the present work. 

An important requirement when performing numerical simulations of black holes in which the singularities are 
excised from the computational domain is to use coordinates that are regular and penetrate the horizon. Foliations 
that penetrate the horizon facilitate the task of removing (i.e. excising) a region containing the black hole singularity 
while preserving the causal structure of the spacetime exterior to the event horizon. Thus the numerical results we 
present here consist of numerical evolutions that for infinite resolution correspond to the solution of a single black 
hole expressed in ingoing-Eddington-Finkelstein (iEF) and also Painleve-GuUstrand (PG) [ p3|j24| . coordinates. The 
iEF coordinates coincide with Kerr-Schild coordinates in the case of zero angular momentum ]25|]. Recently, Martel 
and Poisson |2^ have showed that the PG and iEF coordinates are members of the same one-parameter family of 
coordinate systems. 



In iEF coordinates, the line element (3^) takes the explicit form 



ds^ = - {I ]dt^ + dt dr+ [1 + ]dr^ + r^dfl^ . (3.21) 

r J r \ r 



The ADM variables in these coordinates are given by: 



+ (3.22, 



6 = r (3.25) 
K, = -^{r + m)(l + —) (3.26) 



2m f ^ 2m 



-1/2 

^6 = ^ ( 1 + ^ ) • (3-27) 
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Similarly, for PG coordinates 



and 



ds' = -[l~ ^) + 2 J^dt dr + dr^ + rH^'' , (3.28) 



a = 1 (3.29) 

P = (3.30) 
a = l (3.31) 



6 = r (3.32) 
2r 



i^a = -T^ (3.33) 



Kb = -. (3.34) 



r 



The geometrical interpretation of the iEF coordinate system is that in addition to having a timelike killing vector, 
the combination of timelike and radial tangent vectors dt — dr remains null. In terms of the spacetime metric, this 
condition is stated as gtt — '2gtr + Qrr — 0, or similarly in terms of 3+1 metric functions as a — a{l~ (3). On the other 
hand, the PG coordinate system can be viewed as that anchored to a family of freely moving observers (time-like) 
starting at infinity with vanishing velocity p6t . 

Finally, in numerical evolutions of spherically symmetric spacetimes, it is useful to monitor not only the Hamiltonian 
and momentum constraints but to pay attention also to the mass function 



b 
2 

In vacuum, this M is the gauge invariant definition of mass 



M{r,t)^-{l-V,,bV^b). (3.35) 



IV. ADJUSTED ADM SYSTEMS 



Each of the equations in the g-K or standard ADM system ( |3.9| - p.l2[ ) has the form 

dt W(„) - I3dr -U(„) = S(„) M(„) + C(„) (4.1) 

(no summation over n) with 

U(„) = (a, b, Ka, Kb) (4.2) 

B(„) =dia,g{-aKa + dr/3, -a Kb, 2 a Kb, aKa) (4.3) 

The C(„) contain the remaining terms that cannot be written in the form W(„), with i?(„) independent of U(„). If 



one views each equation as independent, that is fi, and C(„) as given and independent of W(„), Eqs. (4.1) admit 



rapidly growing solutions if i?(n) > 0, exponentially growing if B(n) and j3 are non-negative constants. Obviously, 



when the equations (4.1) are considered as a coupled system of equations, one cannot guarantee the above conclusion. 
However, we have been able to track the problems observed in our numerical evolutions to those terms for which 

B(n) > 0. 

In iEF and PG coordinates, the second and fourth components of i?(n) are non-positive. The first component of 
-B(„) vanishes for PG coordinates and is positive for iEF. Finally, the third component of i?(„) is always non-negative, 
for both iEF and PG coordinates. Our numerical experiments indicate that the origin of the instabilities is due to the 
term involving i3(„=3) — 2 a Kb- The term in iEF coordinates involving i?(„=i) with the "wrong" sign did not seem 
to affect the evolutions (see results below). 

The objective is then to find a way to " change the sign" of i3(„=3). Fortunately, the combination Kb Ka also 



appears in the Hamiltonian constraint ( 3.13| ). One can then add to the evolution equation for Ka a term of the form 



-/z aH. As a consequence, i?(„=3) = 2 {1 — fj,) a Kb- Therefore, in principle, any value of ^ > 1 should produce stable 
evolutions. Our numerical simulations show that = 2 is an optimal value to reach quickly the time independent 
solution. However, the same experiments indicate also that even values of /i « 0.5 yield stable evolutions. The reason 
for this is likely because the simple analysis above does not take into consideration the non-linear coupling in the 
equations. 
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V. LINEAR GAUGE MODES 



A potential source of instabilities in numerical simulations is the presence of gauge modes. Gauge modes arise 
because a prescription of the lapse and shift is not sufficient for complete gauge fixing. If the lapse and shift functions 
are for instance determined from algebraic expressions or differential equations, there is still remaining freedom left 
to perform coordinate transformation that take us from {gij,Kij) to {gij,Kij), leaving the lapse-shift prescription 
invariant. In other words, the data {gij, Kij) will be unique up to coordinate transformations that leave the lapse-shift 
prescription invariant. One can then encounter the situation in which the transformed pair {gij,Kij) continues to 
satisfy the constraints, but it possesses unbounded growths. In principle, because {gij,Kij) are allowed solutions to 
the evolution equations, one could through numerical convergence and monitoring of the constraints single out these 
modes. In practice the situation is not that simple. Gauge conditions that allow the rapid growth of grid functions 
are likely to trigger numerical instabilities. 

To investigate these modes, let us consider an infinitesimal coordinate transformation, 

x'^^x^'-^C^, (5.1) 
This transformation induces a linear perturbation of the spacetimc metric, 

gi^v 9tiy + 5g^_,y , (5.2) 

of the form 

Sgfi,, = -C^9f_i,,- (5.3) 

In spherical symmetry, the most general coordinate transformation is 

t^t + St{r, t) (5.4) 

r ^ r + 5r{r, t) , (5-5) 

or equivalently 

= (5t, (5r, 0, 0) . (5.6) 



Given ( x£ ) , the gauge induced perturbations of the metric takes the form 

5gtt - -2 {-a^ + a^/?^) dtSt - 2a^p dtSr ~ dt {-a^ + a^P^) 6t - dr {-a^ + a^P^) Sr (5.7) 

Sgtr = -a^P dtSt - dfSr - {-a^ + 0^0^) drSt - a^P drSr - dt {a^P) St - dr {a^P) 5r (5.8) 

Sgrr = -2a^P drSt - 2a^ drSr - dt (a^) St - dr (a^) 6r (5.9) 

Sgeo = -dt St - dr (&') Sr (5.10) 

= Sgee sin^e, (5.11) 



with all remaining components vanishing. On the other hand, from ( |3.5| ), the components of the perturbed metric 
Sg^ii, above are also given by: 

Sgtt^ ~2aSa + 2aSaP'^ + 20^ PSP (5.12) 
Sgtr ^2aSaP + a^ SP (5.13) 
Sgrr — '2 a Sa (5-14) 
^2bSb. (5.15) 



In studying these gauge perturbations, we consider a number of different prescriptions for choosing the lapse and 
shift. These impose restrictions on the components of Sgfj_i, to linear order, but are not sufficient for complete gauge 
fixing. 

One piece of evidence usually given to argue for the highly unstable properties of the standard ADM formulation 
is its inability to simulate single black hole spacetimes. Specifically, the problem consists of using a known analytic 
black hole solution of Einstein's equation to construct initial data, set boundary conditions and specify the lapse and 
shift. Given this input, the output is the numerical evolution of the spatial metric gij and extrinsic curvature Kij. 
This simple setup does not yield stable evolutions. 
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Choosing the lapse and shift to be given by the exact analytic solutions implies that in the coordinate trans formation s 



above, Sa = 6(3 = 0. We shall call this the "exact lapse + exact shift" (EL+ES) condition, It is clear from (5.7- 5.11 ) 
that the EL+ES choice does not fix the coordinates uniquely, but rather results in a non-trivial system of PDE's to 
be satisfied by St and Sr. The solutions to this PDE system form an equivalence class of gauges all satisfying the 
same condition on the lapse and shift. At the analytic level, one selects a unique member of this equivalence class by 
imposing boundary and initial conditions. Numerically however, because of truncation errors, this is only possible if 
the evolution is stable 



Substitution oi Sa — 6(3 — in { 5.12 - 5.14 ) yields Sgtr — (3 6grr = 6gtt — (3'^ 6grr — 0. Using ( 5/7 ) , these conditions 



can be rewritten in terms of the gauge perturbation ^'^ as 

dt£, = Adr^ + (5.16) 

where 

^={t) and A. (I I 

The matrix A has a degenerate eigenvalue \ = (3 and corresponding eigenvector e = [0, 1]. Therefore, the system is 
only weakly hyperbolic and thus ill-posed. It is then not possible to guarantee the absence of rapidly growing gauge 
modes. Gauge instabilities, by themselves, do not violate the constraints; however, they are likely to trigger numerical 
instabilities and thereby couple to constraint-violating modes. Note that this conclusion does not depend on the use 
of the standard ADM formulation, but it applies to any initial value formulation of Einstein's equation. This is likely 
the reason why it has not been possible to produce hyperbolic formulations of Einstein equations using the EL-I-ES 
lapse-shift condition. Other groups have reported numerical instabilities associated with an EL+ES prescription. 
(See, for example, |^,^.) Our result, which provides some analytic insight into the source of these instabilities, 
extends the analysis of Ref. |27| by including perturbations of a general spherically symmetric spacetime ( |3.5[) and by 
making no assumptions about the form of the coupling between gauge and constraint-violating perturbations. 

Finally, the inability to guarantee the absence of unbound gauge modes does not necessarily imply that it is 
impossible to design an evolution scheme that is long-term stable and convergent. As we shall see, by adjusting the 
standard ADM system with constraint terms, stable evolutions are possible even in the presence of these gauge modes. 

VI. STABLE ADJUSTED ADM SYSTEMS IN ID 

We consider next a series of lapse-shift choices. For each choice, we investigate (1) the quasi well-posedness of the 
resulting system of adjusted ADM evolution equations, (2) the propagation of the constraints and (3) the convergence 
and stability of numerical simulations. 

A. Exact Lapse + Area Locking 

We can take advantage of the assumed spherical symmetry of the problem and "lock" the area of constant-r surfaces. 



That is, we exploit the lapse-shift freedom and set dt gee — 0, or equivalently dtb = 0\/ t. From Eq. (3.10), this yields 



= -(3drb + abKi, , (6.1) 



which can be seen as an algebraic equation to solve for /3 or a. We will use (6.1) as 



with 6 determined by the initial data. In our case, for both iEF and PG coordinates, b = r. In addition, we choose an 
exact lapse, i.e. an arbitrary but a priori specified function of spacetime. Here again, since the goal is to reproduce 
numerically the analytic solution, we set the lapse to that given by the iEF or PG solutions. 

This exact lapse, area locking (EL+AL) gauge condition has previously been investigated in Ref. p7||. However, 



the implementation of locking the area was done at the numerical level. That is, condition (6.2) was not explicitly 
used. Instead, during the temporal updating of grid functions, a correction to the shift was introduced to keep the 
area locked. With this numerical area-locking and with blending outer boundary conditions, stable simulations were 
reported in Ref.. p7[| for computational domains with ro < 11 to. 
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Interestingly enough, this EL+AL choice of lapse and shift yields an ADM system of equations already first order 
in space and quasi-linear; namely, 



dtu = Adr 



B'i 



(6.3) 



with 




and A 



/3 

r dr a + 2 {1 — fi) a 



V 



aar\ 
P 
(3 J 



where (3 is given by (6.2) and we have set & = r to simplify notation. In the numerical evolutions, however, we do not 
set b ^ r. The numerical code includes the evolution equation for b. We do not explicitly write the matrix B since 
it is not necessary for the analysis below. However, we must emphasize that, in order to achieve stability, a value of 
> 1 is in principle needed (see Sec. IV). 



Although this EL+AL case is not representative of the general ADM equations, where second spatial derivatives 
do appear, we will use it to introduce the main techniques and ideas regarding quasi well-posedness. We notice first 
that the matrix A has eigenvalues 



Ai = /3 

A2 = /? + - 
a 

a 



and corresponding eigenvectors 



ei 



62 



63 



[0,1,0] 



0^ r. 



r 9r a + 2 (1 — /i) a 



2 r dr a + 2 (1 — n) a 
-a r, ,1 



(6.4) 
(6.5) 

(6.6) 

(6.7) 
(6.8) 

(6.9) 



Because all of the eigenvalues are real and distinct, the system ( |6.3| ) is strictly hyperbolic independent of the addition of 
the constraint term. The EL+AL system of equations is then an example of a hyperbolic system that, unless suitable 
constraint-terms are added, is subject to developing rapidly growing solutions (see Ref. B for another example). It 
is also important to notice that Ai represents a characteristic speed corresponding to propagation along the timelike 
normal to the foliation. Similarly, A2 and A3 are characteristic speeds along the light cone. 

Let's now consider the particular case of initial data and lapse function constructed from the single black hole 
solution in iEF or PG coordinates. The eigenvalues are given in these coordinates by 



Ai 
A2 
A3 



2m 



_ I2m 
r + 2m V r 



1 



1 



2 m 
r 

2171 ~ r 1 2 m 



2 m + 7' 



1, 



(6.10) 
(6.11) 
(6.12) 



where the first and second equalities are for iEF and PG coordinates, respectively. Since the excision boundary is 
by construction inside the black hole ho rizo n (i.e. Ve < = 2 m), we have that all the eigenvalues are positive there. 
By looking at the principal part of Eq. (6.3), it is easy to see that non-negative eigenvalues imply propagation of field 
variables in the direction of decreasing r-coordinate. Therefore, at Tg all the fields propagate out of the computational 
domain into the hole singularity, thus boundary conditions are not required. On the other hand, at the outer boundary 
To, one has A3 < 0, with the remaining eigenvalues still non-negative. At r,,, one has then two ingoing modes (Ai and 
A2) and only one outgoing (A3). Given this information, it is perfectly possible to impose a condition suppressing 
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modes entering the computational domain. However, as we mentioned before, we choose not to do so and impose at 
To the analytic, exact solutions. The direct consequence will be the appearance of a pulse at the outer boundary of 
the computational domain due to discontinuities in the truncation errors between the outermost evolved point and 
the boundary point Tq. This pulse propagates in the direction of the black hole and leaves the computational domain 
through the excision boundary. 

The next step is to analyze the effect that the EL+AL choice has on the propagation of the constraints. At the 
continuum level, for arbit rary choice of lapse, shift and initial data (g^j, Kij) satisfying the constraints ( ^.3| ) and 
the evolution equations (3.1) and ( |3.2[ ) guarantee, ignoring for the moment boundary conditions, that the evolved 
data will continue to satisfy the constraints. If one now takes into consideration boundary conditions, it is important 
to keep in mind that boundary data {gij^ Kij) must satisfy the constraints. By looking at the way the constraints 
propagate, i.e. their characteristics, one gains insight about the allowed boundary conditions consistent with the 
constraints. Another important aspect of well-posedness in the propagation of the constraints is that it guarantees 
that there will be no unbounded high frequency growth appearing in the constraints if they are not exactly satisfied 
at the initial slice (for example, due to numerical errors). This well-posedness for the constraint propagation is a 
non-trivial property, not possible to prove for generic formulation of Einstein's equations . 



What we look for are evolution equations for the constraints, equations that would hold if the system (6.3) is 
satisfied. They can be found by taking time derivatives in both s ides of equations ( 3.13 ) and ( |3.14 ), replacing the 
time derivatives of the metric by the r.h.s. of equations (3.9-3.12), and finally expressing the metric and its spatial 
derivatives in terms of the constraints and their spatial derivatives. Following this procedure, it is not too difficult to 
show that 



dtv — P drV 



(6.13) 



where now 



H 

M 



and P 



f 












a 




[ 


4 



(6.14) 



The matrix P has eigenvalues Ai = A2 , A2 = A3 , with A2 and A3 given by (6.5) and ( p.6| ), respectively. This implies 



that the system (3.13) is also strictly hyperbolic with characteristic speeds along the light cone. Notice that the 
constraints at r^, propagate out of the computational domain into the hole singularity, consistent with the outgoing 
propagation of the field variables at re, namely the tilting, into the black hole, of the light cone. At the outer boundary 
To, there is an ingoing mode (A2 > 0). Therefore, one has, as expected, to be careful to provide data at Tq consistent 
with this entering mode. Since we are imposing at the outer boundary the analytic iEF and PG solutions, the data 
at To already satisfy the constraints. However, as mentioned before, choosing the lapse and shift does not completely 
fix the gauge freedom, thus one still has to be careful handling the gauge modes described in Sec. 0. 

Figure ^ shows the L2 no rm of the Hamiltonian constraint as a function of time. The initial data is given by the 
iEF analytic solution (3.21), and the lapse and shift are chosen from the EL+AL conditions. Similar results were 
obtained with PG coordinates. The computational domain extends from r^. ~ Im to Tq — 40 m. We tried larger and 
smaller values for Tq- However, the stability of the simulations was not affected by the location of the outer boundary. 
We use an upwind parameter q — 0.5 and a constraint-term parameter fi = 2. We show runs for resolutions of 
Ar — m/5, to/10, to/20, to/40 with At — 0.25 Ar. The run with Ar = to/5 has a resolution similar to those used in 

three-dimensional simulations of black hole collisions. 

Figure ^ shows the L2 norm of the Hamiltonian constraint (top) and the L2 norm of the mass (3.35) error (bottom) 
taken at time t = 200 to as a function of resolution. The convergence rate implied by the Hamiltonian constraint is 2.18 
and by the mass error 1.7. The convergence rate from the Hamiltonian constraint is larger than second order because 
we used third order accurate discretizations of the advection term as well as third order accurate extrapolations at the 
excision. On the other hand, the convergence rate obtained from the mass error is less than second order because the 
mass function is proportional to g^^, a quantity difficult to handle numerically near the singularity. The reason for 
using the Hamiltonian constraint and mass to monitor accuracies and convergence is because of the gauge invariant 
nature of these quantities. Finally for reference, we show in Fig. ^ the L2 norm of the solution Hamiltonian constraint 
as a function of time for runs with resolution Ar = to/10 varying the parameter fi. Different lines correspond to 
values of /i = 0.0, 0.025, 0.5, and 2 in order of stability improvement. It is clear from this figure the dramatic effect 
that the added constraint term has on the stability of the simulations. 
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B. Ingoing Null + Area Locking 



The ingoing null + area locking (IN+AL) recipe to specify the lapse and shift consists of, in addition to locking 
the areal coordinate, imposing the condition that the vector dt ~ dr remains null throughout the evolution. This 
null condition is only compatible with the iEF case since by construction the iEF solution is based on ingoing null 
observers. An analogous (ingoing timelike) condition can be obtained for the PG case. In terms of the spacetime 
metr ic, the ingoing null condition is stated as gu — gtr + .9rr = 0, or similarly in terms of 3+1 metric functions in 
(3.5) as 



Conditions ( 6.15| ) and (6.1) yield 



arKb + 1 



a (1-/3). 



and (3 — 



arKh 



arKb + 1 



3.15) 



(6.16) 



where we have set b = r since by construction b remains locked to r. This prescription for the lapse and shift has 
been successfully applied in the past pO|-^. 

In order to make the ADM equations in the IN+AL gauge a first order in space system, we need to introduce 
two new variables: w = dr a and y = dr Kb- Even after this choice is made, there is no unique way of writing the 
resulting equations as a quasi-linear system. The reason is the ambiguity one encounters when dealing with the terms 
involving dr a. One has the choice to either keep it as dr a, substitute it with w, or a combination of both. Either 
choice changes the principal part of the equation. However, in our case, we need only to find a choice that yields a 
strongly hyperbolic system. It turns out that the simplest choice of replacing dr a hy w everywhere yields a system 
that is well-posed. That is, the resulting ADM equations have the form 



with 



dtu = AdrU + Bu, 



( \ 



(6.17) 



Kb 

y 

V w ) 



(6.18) 



and 



/ 

zar Kb 




zaKb z + aV^if? 





a-^Kb 



(6.19) 



VO -za^ c? r ar Kb{l + z) J 

where we have introduced z = I + ar Kb to simplify notation. The corresponding eigenvalues and eigenvectors are 

(6.20) 
(6.21) 
(6.22) 
(6.23) 

(6.24) 



Ai = with eigenvectors ei = [1, 0, 0, 0, 0] , 62 = [0, 0, 1, 0, 0] , 
X2 — 1 with eigenvectors 63 — [0, 1, 0, 0, — za^] , 64 — [0, 0, 0, 1, a^r] 
arKb — 1 



A, 



with eigenvector 65 = [0, 1, 0, ^aKb, a ] 



Notice that all the eigenvectors are independent, and, thus, the system is strongly hyperbolic. Also notice that 
the eigenvalues A2 and A3 are again the characteristic speeds {f3 ± a/a) along the light cone. Furthermore, in iEF 
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coordinates, A3 = (2 m — r)/(2m + r). Therefore, one encounters a similar situation to that of the EL+AL case; 
namely, at re all the eigenvalues are non-negative, and at Tq one has that only A3 < 0. The existence, uniqueness, 
and well-posedness for IN+AL follows then as with EL+AL case. 

Regarding the constraints, their evolution is also described by a strongly hyperbolic system with characteristic 
speeds along the light cone. It is important to stress that our ADM equations already imply this, i.e. it does not 
have any relation to making the evolution equations first order in space. In fact, the principal part of the evolution 
equat ions for the constraints is exactly the same as in the EL+AL gauge, but now the lapse and shift are given by 
( 5.16 ). Thus, the analysis and conclusions also follow as in that case. 

An interesting aspect of the IN+AL choice is that it is possible to find the general solution to Einstein's equations. 
We start by defining / = raKt,. We then use the momentum constraint ( |3.14 ) to e li minate Kg from the other 
equations. The outcome is that we nee d on ly to solve t hree of th e fou r equations ( |3.9| ), ( p.ll| - pTl3 ). We choose to 
work with the Hamiltonian constraint, (3.13), and Eqs. (3.9) and (3.12). The resulting system of equations reads 



0^a{f + -l) + 2rafdrf + 2r{l- f)dra 



{dt - dr) a 



{dt - dr) f 



1) 



2r(l + /)2 



^(1 + /) 



The general solution of these equations is 



2 ra 
r 



C I 



2 m 
r 



a' = (l + C)(l + /) 

where C ^ C{t + r) and m a constant. In terms of the 3+1 variables, the solutions reads 



a2 = (l + C)(l + /) , b^^r^ , Ka=dr 



f 



Kb^ — , a 

ar 1 + / 



with the line element (3.5) explicitly given by 

2 m 



(l + Cy [ I ~ — ] dt^ + 2(1 + C) 



2 m 
r 



1-^)C 
r 



dt dr 



1 + /' 



2 m 



1-^|C 
r 



dr^ + dn^ . 



(6.25) 
(6.26) 

(6.27) 



(6.28) 
(6.29) 

(6.30) 

(6.31) 
(6.32) 



By setting C = one r ecovers the iEF solution ( 3.21 ). Also, it is not difficult to show from the gauge invariant 
definition of mas s (|3.35| ) that the parameter m is indeed the mass of the black hole. An important property of the 
general solution ( |6.31| ) is that it shows explicitly the residual gauge freedom associated with the IN+AL choice of 
lapse and shift. We have found explicitly the equivalence-class of solutions gauge related to the iEF solutions that 
satisfy the IN+AL lapse-shift condition. 

We repeat the same type of numerical experiment as with the EL+AL case, same parameter val ues, r esolutions, 
boundary conditions and initial data. The lapse and shift are however constructed in this case from ( |6.16|) . Figure ^ 
shows the L2 norm of the Hamiltonian constraint as a function of time for different resolutions. Figure H shows the 
L2 norms of the Hamiltonian constraint (top) and mass error (bottom) for different resolutions. The convergence rate 
are similar to those in EL+AL, namely 2.19 from the Hamiltonian constraint and 1.7 from the mass error. Finally, 
Fig. ^ shows the L2 norm of the Hamiltonian constraint as a function of time for a resolution of Ar = m/10. Different 
lines correspond to values of ^ = 0.0, 0.025, 0.5, and 2 in order of stability improvement. It is clear from this results 
that the stability behavior of the system of equations under the IN+AL gauge choice follows closely that of EL+AL. 



C. Exact Lapse + Exact Shift 

Finally, we consider the case in which the lapse function and shift vector are prescribed by the exact analytic 
solutions. To date, it has not been possible to obtain with the standard ADM system long-term stable and convergent 
numerical evolutions under the EL+ES choice. Before we present results of EL+ES evolutions using the adjusted 
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ADM system of equations, let us investigate its quasi well-posedness properties. To make the evolution equations first 
order in space, we add a new variable y = drb. Once more, the resulting set of equations have the form 



dtu^AdrU-^-Bu, 



(6.33) 



with 



and 



I " \ 

6 

Kb 

V y J 



P 





(6.34) 



b drOi + 2 (1 — /i) ay 







ay 





The matrix A has eigenvalues 



Ai 

A2 

As 
A4 
As 



p 

0/3 
-Q 



P 


P + a/a 
P — a/a . 






2 (1 - ^)a 
a 



(6.35) 



(6.36) 
(6.37) 
(6.38) 
(6.39) 
(6.40) 



There are two fields that propagate with characteristic speed (P) along the normal to the hypersurfaces in the foliation, 
one field with zero speed and two other fields with characteristic speeds {P ± a /a) along the light cone, ft is not 
difficult to show that the eigenvectors corresponding to the eigenvalue Xi = X2 — P are not distinct, whatever the 
value of /X. Therefore, the system of equations is only weakly hyperbolic and thus not quasi well-posed, as we have 
already seen by considering the EL+ES gauge perturbations by themselves in Sec. 0. 
As in the previous two cases, the the constraints propagate according to 



where now 



and 



dtV — P drV + I 



H 
M 



(6.41) 



(6.42) 



P 



This is, we have the same situation as before, namely propagation of the constraints along with characteristic speeds 
along the light cone. 

We have carried out the same numerical experiments as for the previous two cases. However, for the iEF solution, 
the simulations with computational domains for Tq > 40 to crashed. We have been able to identify two possible 
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sources behind this problem. One of them is the pulse originated at the outer boundary due to discontinuities of 
the truncation error. This pulse propagates in the direction of the black hole (i.e. decreasing r-coordinate) with the 
characteristic speed [3. As the pulse travels its amplitude grows in time. This effect is shown in Fig. |^. Here we plot 
the solution error for the metric function a as a function of space and time for a resolution Ar = m/10 and outer 
boundary located at — 40 m. The simulation stops because this pulse-error increases to the point that the metric 
fmiction a becomes negative. For small computational domains [tq < 40 to), the crossing time of this pulse is short 
enough and does not allow a catastrophic growth of the pulse. The pulse is able to leave the computational domain 
through the excision boundary without crashing the simulation. Since the initial amplitude of this pulse is O(Ar^), 
i.e. the accuracy of the discretization, in principle one could find a resolution small enough such that the growth 
of the pulse does not affect the life of the simulation. However, accessing those fine resolutions in three-dimensional 
simulations is likely to be impractical. The second, and perhaps more severe, source of the problem is the pr esenc e 
of the zero velocity mode. Zero velocity modes were in principle also allowed in the case of IN+AL, see Eq. (3.20). 
However, the adjusted IN-I-AL system yields only numerical solutions of the general form (6.31), which clearly does 
not contain a zero velocity mode. If /i = 0, these zero velocity modes are not suppressed and eventually terminate the 
simulation. This catastrophic effect induced by zero velocity modes has been previously noticed by Alcubierre et.al 

Figure g| shows the L2 norm of the Hamiltonian constraint as a function of time with resolution Ar — m/10 
for different locations of the outer boundary. For computational domains with approximately To < 40 to, the zero 
velocity mode is still present, but is damped eventually. The reasons why this mode stops growing remain unclear. 
Nonetheless, there is a strong indication that this behavior is connected to the particular choice of coordinates used 
to set the exact lapse and shift. If one, instead of the iEF solution, sets the lapse and shift from the PG solution, the 
outcome of the simulations is completely different. EL-I-ES simulations with PG lapse and shift are long term stable 
and convergent for arbitrary sizes of the computational domain as long as the system of equations was adjusted with 
/X > 1. Figure ^ shows the L2 norm of the solution error for the metric function a as a function of time for different 
resolutions with Tq = 40 to for EL-I-ES in PG coordinates. 



VII. CONCLUSIONS 



We have demonstrated that, at least for the case of single black hole spacetimes in spherical symmetry, it is possible 
to obtain long-term stable and convergent numerical simulations using the standard ADM system of equations if 
the equations are adjusted by introducing terms involving the constraints. Results were presented for three choices 
of lapse and shift. In addition, we introduced the concept of quasi well-posedness, which appears to be useful in 
characterizing the properties of the system of evolution equations. We are currently investigating the extension of 
this approach to three-dimensional evolutions. 
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FIG. 1. Z/2 norm of the Hamiltonian constraint as a function of time for iEF initial data and EL+AL lapse-shift. 
The computational domain extends from re = 1 m to ro = 40 m. Lines from top to bottom correspond to resolutions of 
Ar = m/5, m/10, m/20, m/40, respectively. 
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FIG. 2. L2 norm of the Hamiltonian constraint (top) and L2 norm of the mass error (bottom) as a function of resolution. 
The errors plotted were obtained at t = 200 m. The convergence rate implied by the Hamiltonian errors is 2.18 and by the 
mass error 1.7. 
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FIG. 3. Z/2 norm of the Hamiltonian constraint as a function of time for A = m/10 and To = 40 m. Each line corresponds to 
different values of the constraint-term parameter p. The values of p are 0.0, 0.025, 0.5 and 2 in order of stability improvement. 




FIG. 4. Same as in Fig.m but for the IN+AL case. 
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FIG. 5. Same as in Fig.| but for the IN+AL case. 
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FIG. 6. Same as in Fig.| but for the IN+AL case. 



18 



FIG. 7. Solution error for the metric function a in iEF coordinates for the EL+ES case. The resolution is Ar = m/lO, and 
the outer boundary is located at Vo = 40 m. A pulse originated at the location of the outer boundary due to discontinuities in 
the truncation error propagates in the direction of the black hole (decreasing r-coordinate). The simulation stops because this 
pulse-error increases to the point that the metric function a becomes negative. 
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FIG. 8. L2 norm of the Hamiltonian constraint in iEF coordinates for the EL+ES case as a function of time. 
The resolution is Ar = m/10 and fj, = 2. Each line corresponds to different location of the outer boundary: 
Vo = 20m (solid), 30 m (short dash), 40 m (long dash), 50m (short dash long dash). 
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FIG. 9. 1/2 norm of the solution error in PG coordinates of the metric function a as a function of time. Each line correspond 
to different resolutions, Ar = m/5, m/10, m/20, m/40, respectively from top to bottom. 
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